/*!  \file  TriEigen.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      07/02/2013 04:13:06 PM
 *   \copyright GNU Public License.
 */
#include "TriEigen.h"
#include "Matrix.h"
#include <stdio.h>

TriEigenUnitTest::TriEigenUnitTest(TriEigen & ev_): ev(ev_) {
	int l = strlen(ev.getName());
	char buf[l + 20]; 
	sprintf(buf, "TriEigen %s", ev.getName());
	setName 				(buf);
	setDescription  	("Test symetric tridiagonal eigen value implementation");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	setCompileTime(__DATE__, __TIME__); // each class that inherits this one should also have this line.
	setCompileTime(ev.getCompileTime());
}
bool TriEigenUnitTest::checkNorms( const Matrix & T, const Matrix & S, const Matrix & Q, const double limit) const {
	bool res = true;
	const int  nNorms = 5;
	double norm[nNorms];
	norm[0] = (T * Q - Q * S).nrm2();
	norm[1] = Q.rows() <= Q.cols() ? (T - Q * S * Q.T()).nrm2() : 0;
	norm[2] = (Q.T() * T - S * Q.T()).nrm2();
	norm[3] = Q.rows() >= Q.cols() ? (Q.T() * T * Q - S ).nrm2() : 0;
	norm[4] = Q.OrthNorm(); 
	for(int i = 0; i < nNorms; i++) {
		res &= norm[i] < limit;
	}
	return res;
}
void TriEigen::computeRange( const Matrix & T, Matrix & Sigma, Matrix & Q, const int indexLow, const int indexUp) {
	MISSING;
}
bool TriEigenUnitTest::test() const {
	const unsigned int n = 10;
	double d [] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1.0/8, 2.0/8, 3.0/8,
						4.0/8, 5.0/8, 6.0/8, 7.0/8, 8.0/8, 9.0/8};
	double s_ref [] = {   9.841229182759271e-01,   
	1.952368754827842e+00,   2.920614591444129e+00,   3.888860451192102e+00,
   4.857109715967714e+00,   5.825575957905979e+00,   6.799771278994464e+00,
   7.830775258011721e+00,   9.084898199079339e+00,   1.085590287430078e+01};
	Matrix T(n, n, TriDiagonal, d);
   Matrix S(n, n, Diagonal,d);
   Matrix Sigma_ref(n, n, Diagonal, s_ref);
	Matrix Q(n, n, Orthogonal);	
	ev.compute(T, S, Q);
	bool res = ev.computaionSuccessed();
	const double limit = 1e-13;
	res &= checkNorms(T, S, Q, limit);
	
	T.setView(2, 2, 5, 5);
	Q.setView(0, 0, 5, 5);
	S.setView(3, 3, 5, 5);

	ev.compute(T, S, Q);
	res &= checkNorms(T, S, Q, limit);

	T.setView(2, 2, 5, 5);
	Q.setView(0, 0, 5, 2);
	S.setView(3, 3, 2, 2);

	if (ev.canComputeRange) { // test it :) 
		ev.computeRange(T, S, Q, 1, 2);
		res &= checkNorms(T, S, Q, limit);
	}

	return res;
}
